/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

#include "definitions.h"

// Standard Lennard-Jones potential 
void CalcLennJonesInteraction( real *ijEnergy, // out - interaction energy
			       VecR *ijForce, // out - force vector
			       const real ijSigSq, 
			       const real ijEps,
			       const real rr,
			       const VecR r )
{
  real sr2, sr6, sr12;
  
  sr2 = ijSigSq / rr;
  sr6 = Cube( sr2 );
  sr12 = Sqr( sr6 );
  *ijEnergy = 4. * ijEps * ( sr12 - sr6 ); 
  VSCopy( *ijForce, 24. * ijEps * ( 2. * sr12 - sr6 ) / rr, r );
}

// Shifted-Force LJ potential [Stoddard and Ford, Phys Rev A, 8, 1504 (1973)]
void CalcLennJonesInteraction_SF( real *ijEnergy, // out - interaction energy
				  VecR *ijForce, // out - force vector
				  const real ijSigSq, 
				  const real ijEps,
				  const real rr,
				  const real rrCut,
				  const VecR rij )
{
  real sr2, sr6, sr12, srC2, srC6, srC12, t;
  
  sr2 = ijSigSq / rr;
  sr6 = Cube( sr2 );
  sr12 = Sqr( sr6 );
  srC2 = ijSigSq / rrCut;
  srC6 = Cube( srC2 );
  srC12 = Sqr( srC6 );
  t = 2*srC12-srC6;
  *ijEnergy = 4 * ijEps * ( sr12 - sr6 + 3*t*srC2 - 7*srC12 + 4*srC6 );
  VSCopy( *ijForce, 24 * ijEps * ( ( 2*sr12 - sr6 ) / rr - t / rrCut ), rij );
}

